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to study the numerical evolution of D dimensional vacuum space-times with an SO(D — 2) isometry 
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3+1 dimensional numerical codes with small adaptations. Brill-Lindquist initial data are constructed 
in D dimensions and a procedure to match them to our 3+1 dimensional evolution equations is 
given. We have implemented our framework by adapting the Lean code and perform a variety 
of simulations of non-spinning black hole space-times. Specifically, we present a modified moving 
puncture gauge which facilitates long term stable simulations in D = 5. We further demonstrate the 
internal consistency of the code by studying convergence and comparing numerical versus analytic 
results in the case of geodesic slicing for D — 5, 6. 
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1. INTRODUCTION 



Numerical relativity is an essential tool to study many processes involving strong gravitational 
fields. In four space-time dimensions, processes of this sort, such as black hole (BH) binary 
evolutions, are of utmost importance for understanding the main sources of gravitational waves, 
which are expected to be detected by the next generation of ground based [Laser Interferometer 
Gravitational- Wave Observatory (LIGO), VIRGO] and space based [Laser Interferometer Space 
Antenna (LISA)] interferometers. Long-term stable numerical evolutions of BH binaries have 
finally been achieved after four decades of efforts [H-Q. The numerical modelling of generic 
spinning BH binaries in vacuum Einstein gravity is an active field of research, with important 
consequences for gravitational wave detection in the near future. 

Numerical relativity in a higher dimensional space-time, instead, is an essentially unexplored 
field, with tremendous potential to provide answers to some of the most fundamental questions 
in physics. Recent developments in experimental and theoretical physics make this a pressing 
issue. We refer, in particular, to the prominent role of BHs in the gauge-gravity duality, in 
TeV-scale gravity or even on their own as solutions of the field equations. These are some of the 
most active areas of current research in gravitational and high energy physics. 



l.l. 



Motivation 



i. AdS/CFT and holography. In 1997-98, a powerful new technique known as the AdS/CFT 
correspondence or, more generally, the gauge-string duality, was introduced and rapidly 
developed Q. This holographic correspondence provides an effective description of a 
non-pcrturbativc, strongly coupled regime of certain gauge theories in terms of higher- 
dimensional classical gravity. In particular, equilibrium and non-equilibrium properties of 
strongly coupled thermal gauge theories are related to the physics of higher-dimensional 
BHs, black branes and their fluctuations. These studies revealed intriguing connections 
between the dynamics of BH horizons and hydrodynamics Q , and offer new perspectives 
on notoriously difficult problems, such as the BH information loss paradox, the nature of 
BH singularities or quantum gravity. 

Numerical relativity in anti-dc Sitter backgrounds is bound to contribute enormously to 
our understanding of the gauge-gravity duality and is likely to have important applications 
in the interpretation of observations 6 9] . For instance, in the context of the gauge-gravity 
duality, high energy collisions of BHs have a dual description in terms of a) high energy 
collisions with balls of de-conhncd plasma surrounded by a confining phase and b ) the rapid 
localised heating of a de-confined plasma. These are the type of events that may have 
direct observational consequences for the experiments at Brookhaven's Rclativistic Heavy 
Ion Collider (RHIC) HQ. Numerical relativity in anti-de Sitter is notoriously difficult, 
and so far only very special situations have been handled [l(| [H| • The phenomenologically 
most interesting case is a five dimensional space-time, AdS§, and therefore the higher 
dimensional extension of numerical relativity is necessary. 

ii. TeV-scale gravity scenarios. An outstanding problem in high energy physics is the ex- 
tremely large ratio between the four dimensional Planck scale, 10 19 GeV, and the elec- 
troweak scale, 10 2 GeV. It has been proposed that this hierarchy problem can be resolved 
if one adopts the idea that the Standard Model is confined to a brane in a higher dimen- 
sional space, such that the extra dimensions are much larger than the four dimensional 
Planck scale (they may be large up to a sub-millimetre scale) • I n a different version 

of the model, the extra dimensions are infinite, but the metric has an exponential factor 
introducing a finite length scale [HI, [l6| . 

In such models, the fundamental Planck scale could be as low as 1 TeV. Thus, high energy 
colliders, such as the Large Hadron Collider (LHC), may directly probe strongly coupled 
gravitational physics [17H22j |. In fact, such tests may even be routinely available in the 
collisions of ultra-high energy cosmic rays with the Earth's atmosphere [23l425| , or in astro- 



thc production of BHs at trans-Planckian collision energies (compared to the fundamen- 
tal Planck scale) should be well described by using classical general relativity extended 
to D dimensions [T^-HM I29h33ij. The challenge is then to use the classical framework to 
determine the cross section for production and, for each initial setup, the fractions of the 
collision energy and angular momentum that are lost in the higher dimensional space by 
emission of gravitational waves. This information will be of paramount importance to im- 
prove the modelling of microscopic BH production in event generators such as Truenoir, 
Charybdis2, Catfish or Blackmax pol l33 - [37| . The event generators will then provide 
a description of the corresponding evaporation phase, which might be observed during LHC 
collisions. 

The first models for BH production in parton-parton collisions used a simple black disk 
approach to estimate the cross section for production [llj [2(J. Improved bounds have 
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been obtained using either trapped surface methods to estimate the cross section for BH 
production [38l - l4l| or approximation schemes (42W471 ] to evaluate the gravitational energy 
loss. Only recently exact results for highly rclativistic collisions where obtained in four 
dimensions, using numerical relativity techniques j48l - l50j . No such exact results are yet 
available in the higher dimensional case. To obtain them is one of our main goals and the 
present paper introduces a formalism to achieve that. 

Higher dimensional black holes. Asymptotically flat higher dimensional black objects have 
a much richer structure than their four dimensional counterparts. For instance, spherical 
topology is not the only allowed topology for objects with a horizon. One can also have, e.g., 
black rings, with a donut-like topology. Remarkably, these two different horizon topologies 
coexist for certain regions in phase-space [5l| . The stability of general higher-dimensional 
BHs is now starting to be explored. Generically it has been conjectured that for D > 6 
ultra-spinning Myers- Perry BHs will be unstable [52| . This instability has been confirmed 
by an analysis of linearised axi-symmetric perturbations in D = 7, 8, 9 [53j . Clearly, the 
study of the non-linear development of these instabilities requires numerical methods, such 
as the ones presented herein. A study of this type was very recently presented for a non 



axi-symmetric perturbation in D = 5 [541 ]. where it was found that a single spinning five 
dimensional Myers-Perry BH is unstable, for sufficiently large rotation parameter (thereby 
confirming previous conjectures [55rl57^ . 

Not much is known about general equilibrium states in anti-de Sitter backgrounds. The 
gauge-gravity duality and the hydrodynamic limit have been used to predict the existence of 
larg er classes of BHs in anti-dc Sitter backgrounds, including non axi-symmetric solutions 
[56l \5l\. However, these have not yet been found. 

Finally, there are issues of principle, as for example testing Cosmic Censorship in BH collisions 
pel l5(i| which require state-of-the-art numerical simulations. 



1.2. Space-times with symmetries 



From what has been said, the extension of four dimensional numerical Relativity is mandatory. 
Some pioneering works have been concerned with the non-linear development of the Gregory- 
Laflamme instability [f58| of cosmic strings [59| and gravitational collapse, with spherical symme- 
try jr30| , axial symmetry [6lT ] or even static situations [62| . Another numerical code, based on the 
cartoon method [63[, was developed and tested for five space-time dimensions in Ref. [64| . See 
also Ref. [(35| for a discussion of slicings of D dimensional black holes. The (phenomenologically) 
most interesting large extra dimensions models are, however, in higher than five space-time di- 
mensions (see for instance [30|V Moreover, the ultra-spinning instabilities of Myers-Perry BHs 
should occur in D > 6. Thus, our approach here is to develop a framework and a numerical code 
that can, in principle, be applied to different space-time dimensions with little adaptations. This 
may be achieved by taking the D dimensional vacuum space-time to have an isometry group 
fit to include a large class of interesting problems. If this isometry group is sufficiently large, it 
allows a dimensional reduction of the problem to 3+1 dimensions, wherein it appears as (four 
dimensional) general relativity coupled to some quasi-matter terms. 1 Thus, the different space- 
time dimension manifests itself only in the different quasi-matter content of the four dimensional 



Hereafter, we dub the source terms of the lower dimensional Einstein equations as quasi-matter, since its 
energy-momentum tensor is not that of canonical matter. 
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theory. We emphasise, in this context, that full blown 4 + 1, 5 + 1, etc. numerical simulations 
without symmetry are currently not possible due to the computational costs, so that our ap- 
proach pushes numerical relativity in higher dimensions to the outmost practical limits of the 
present time. Moreover, an obvious advantage of this approach is that we can use existing codes 
with small adaptations: the four dimensional equations need to be coupled to the appropriate 
quasi-matter terms and some issues related to the chosen coordinates must be addressed, as we 
shall see. Finally, the lessons learnt in treating our effective gravity plus quasi-matter system 
might be of use in dealing with other four dimensional numerical relativity problems with sources. 

1.3. Axial symmetry SO(D - 2) and SO{D - 3) 

We consider two classes of models, which are generalisations of axial symmetry to higher 
dimensional space-times: a D > 5 dimensional vacuum space-time with an SO(D — 2) isometry 
group, and a D > 6 dimensional vacuum space-time with an SO(D — 3) isometry group. The 
former class allows studies of head-on collisions of non-spinning BHs. In order to end up with a 
3 + 1 dimensional model we use, however, only part of this symmetry: we perform a dimensional 
reduction by isometry on a (D — 4)-sphere which has an SO(D — 3) C SO(D — 2) isometry 
group. The latter class allows to model BH collisions with impact parameter and with spinning 
BHs, as long as all the dynamics take place on a single plane. 2 In this case we perform a 
dimensional reduction by isometry on the entire SO(D — 3) isometry group. This class includes 
the most interesting physical configurations relevant to accelerator — and cosmic ray — physics (in 
the context of TeV-scale gravity) , and to the theoretical properties of higher-dimensional black 
objects (such as stability and phase diagrams). 

We formulate the evolution equations in the Baumgartc, Shapiro, Shibata and Nakamura 
(BSSN) formulation (6(| [67}, together with the moving puncture approach @, Q. This is known 
to provide a stable evolution scheme for vacuum solutions in four dimensions, and therefore it 
is the natural framework for our Einstein plus quasi-matter system. The quasi-matter terms 
however, exhibit a problem for numerical evolution, well known from other numerical studies 
using coordinates adapted to axial symmetry, which is sourced by the existence of a coordinate 
singularity at the axis. In our formulation, this problem appears when a certain 3+1 dimen- 
sional Cartesian coordinate vanishes, y = 0. We present a detailed treatment of this problem, 
introducing first regular variables, then analysing one by one all potentially pathological terms 
in our evolution equations and finally presenting a method to heal all of them. The resulting 
equations have no further (obvious) problems for numerical evolution and could, in principle, be 
implemented in any working 3+1 dimensional numerical relativity code. 

Here we present numerical results using the Lean code [6§| . developed by one of us. We stress 
that the formalism developed here is valid in general D. However, long term stable evolutions 
typically require some experiments with free parameters in the gauge conditions and also possibly 
with constraint damping. For D = 5 we show that, if appropriate gauge conditions are chosen, 
the numerical evolution for Brill-Lindquist initial data describing a single BH is stable and the 
constraints are preserved in the evolution, within numerical error. As another test, we evolve 
the same initial data in a geodesic slicing gauge. This gauge is inappropriate for a long term 
evolution; but it allows us to compare the numerical evolution with the analytic solution for a 
single Tangherlini BH in D = 5. We find excellent agreement between the two. We also present 
some preliminary results for D = 6. 



2 This follows from the fact that the angular momenta of the black holes are parallel to the orbital angular 
momentum. 
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This paper is organised as follows. In Section^ we discuss the D dimensional ansatz, perform 
the dimensional reduction by isometry, perform the Arnowitt-Descr-Misncr (ADM) split and 
present the BSSN formulation of our equations. In Section [3l the construction of Brill-Lindquist 
initial data in D dimensions is discussed and a procedure to match it to our 3+1 formulation 
is given. In Section 2] we present the numerical treatment and results. We draw our conclusions 
and discuss implications of our results for future work in Section [5] A considerable part of the 
technical details for the numerical treatment is organised into three appendices. In Appendix [A] 
we motivate and discuss the introduction of regular variables at y = and present all relevant 
equations in terms of these variables. In Appendix[B]we explain how to tackle all the problematic 
terms at y = in these equations. Finally, in Appendix [Cj we discuss the construction of the 
geodesic slicing which is used to compare analytical with numerical results. 

2. THE EFFECTIVE 3+1 DIMENSIONAL SYSTEM 

The starting point of the formalism used here is a dimensional reduction from D dimensional 
general relativity in vacuum to a four dimensional model. The isometry group of D dimensional 
Minkowski space-time is ISO(l, D— 1); solutions of general relativity (or of other metric theories 
of gravity) generically break this symmetry into a subgroup. For instance, the isometry group 
of a Schwarzschild (or, for D > 4, Tangherlini [69| ) BH is SO(D — 1) x K, whereas for a 
head-on collision of two non-rotating BHs it is SO(D — 2): indeed, neither the time direction 
nor the direction of the collision correspond to symmetries, but a rotation of the remaining 
D — 2 spatial directions leaves the space-time invariant. The total space-time can then be 
considered as the semi-direct product of a three dimensional space-time Af with the sphere 
S D ~ 3 = SO(D — 2)/SO(D — 3). A coordinate system for Af can be given, for example in the 
case of a head-on collision of two BHs, by the time t, the coordinate z along the collision axis, 
and the distance from that axis. 

One can take advantage of this symmetry to reduce the space-time dimensionality. This can 
be accomplished by writing Einstein's equations in D dimensions in a coordinate system which 
makes the symmetry manifest, allowing for a lower dimensional interpretation of the D dimen- 
sional Einstein's equations (in the spirit of Kaluza-Klein reduction). We remark, however, that 
we are not performing a compactification; rather, we perform a dimensional reduction by isom- 
etry, as first proposed by Geroch p70j | . The extra dimensions manifest themselves in the lower 
dimensionality as a source of Einstein's equations, defined on the lower dimensional manifold. 

In principle, one could use the symmetry in a more naive way, assuming that the solution 
does not depend on the coordinates parameterizing the sphere and simply evolving the relevant 
components of the D dimensional Einstein's equations. The perspective provided by dimensional 
reduction, however, has two advantages: (i) all quantities have a geometrical interpretation, and 
this allows for a deeper understanding of the problem and a better control of the equations; (ii) 
it is possible to use, with minor modifications, the numerical codes which have already been 
written to implement Einstein's equations in a four dimensional space-time. Therefore, we do 
not use the entire SO(D — 2) symmetry of the process, but only a SO(D — 3) subgroup. This 
reduces the space-time on a (D — 4)-sphere and yields a four dimensional manifold. 

In the original proposal of Geroch [7(| the symm etry sp ace was SO(2). This approach has 
been applied to numerical relativity, see for instance |7lJ-|73| ; a five dimensional extension, with 
the same symmetry space, has been derived in (74|. A generalisation to coset manifolds (like the 
sphere S n ) was given by Cho in [T5|, [76[, but in these papers the complete form of Einstein's 
equations was not presented. Here we provide the explicit form of Einstein's equations for 
symmetry spaces S n together with their numerical implementation. 
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2.1. 4 + (D -4) split 

We now describe in detail the reduction from D to 4 dimensions. In order to highlight the 
particular classes of BH binaries we are able to study with this framework, it is convenient to 
begin this discussion with the isometry group of the S D ~ 3 sphere, i. e. with the 3 + (D — 3) split. 

A general D dimensional space-time metric may be written in the form 

ds 2 = g MN dx M dx N = g-^{x M )dx^dx v + %(x M ) (dx 1 - A l -^x M )dx^ (dx 3 - Ai(x M )dx^ , 

(2.1) 

where we have split the space-time coordinates as x M = (x^,x % ); M,N = 0, ...,D— 1 are 
space-time indices, /2, v = 0, 1, 2 are three dimensional indices and i, j = 3, . . . D — 1 are indices 
in the remaining D — 3 dimensions. We may think of the space-time as a fibre bundle; {x 1 } arc 
coordinates along the fibre and {x^} are coordinates on the base space. 

We are interested in studying D dimensional space-times with an SO(D — 2) isometry group. 
This is the isometry group of the S D ~ 3 sphere, which justifies why we are performing a 3 + (I? — 3) 
splitting of the D dimensional space-time. Thus, we assume that £ a , a = 1, . . . , (D — 3)(D — 2)/2, 
are Killing vector fields, 

^ a 9MN = , (2.2) 

with Lie algebra 

=ea6 C £ c , (2.3) 

where e a b c are the structure constants of SO(D — 2). Because the fibre has the minimal dimension 
necessary to accommodate (D — 3)(D — 2)/2 independent Killing vector fields, we may assume 
without loss of generality that the Killing vector fields have components exclusively along the 
fibre: £ a — t^a&i. Furthermore, we may normalise the Killing vectors so that they only depend 
on the coordinates of the fibre, i.e. d^l = 0. Then Eq. (|2.2[) gives the following conditions 

A- <-->,.. (2-4) 

^4 = °> ( 2 - 5 ) 

Cu9fiP = 0- (2-6) 

These expressions can be interpreted either as Lie derivatives of rank-2 tensors defined on the 
D dimensional space-time, or as Lie derivatives of a rank-2 tensor, a vector and a scalar, which 
are defined on S D ~ 3 . 

Conditions (|2.4p - (|2.6[) have the following implications: 



because, from (|2.4p . f^j admits the maximal number of Killing vector fields and thus must 
be the metric on a maximally symmetric space at each x* 1 . Due to (|2.3p this space must 
be the S D ~ 3 sphere, h?-, denotes the metric on an S D ~ 3 with unit radius; 



9Hu = QaM) , (2.8) 

because the Killing vector fields £ a act transitively on the fibre and therefore the base space 
metric must be independent of the fibre coordinates; 



8 



4=0, (2.9) 

because Eq. (12. 5[) is equivalent to 

Ka,^]=0, (2.10) 

and there exist no non-trivial vector fields on S D ~ 3 for D > 5 that commute with all Killing 
vector fields on the sphere. 

We remark that (|2.10p corresponds to the statement, expressed in [75[ in group theoretical 
language, that the gauge group for a theory reduced on a coset space G/H is the normaliser of 
H in G; in the case of a sphere, where G = SO(D — 2) and H = SO(D — 3), the normaliser 
vanishes and then there are no "gauge vectors", i.e., no non- vanishing metric components g^. If 
the normaliser of H in G is non- vanishing, such metric components appear, and with dimensional 
reduction they yield vector fields which contribute to the stress-energy tensor in the reduced 
theory. For example, in the case of head-on collision, if D = 4, the isomctry space is SO(2) and 
the quasi-matter of the reduced theory consists of a scalar field and of a vector field (as in [70l | 
and in [7ll473] ]): if D > 4, the isometry space is SO(D — 2)/SO(D — 3), and the quasi-matter 
of the reduced theory consists of a single scalar field. In the remainder of this work we focus on 
this subclass of space-times, which already contains a vast class of physically relevant problems, 
and postpone a discussion of the general case with A 1 ^ ^ (i.e., with g^ ^ 0) to future work. 

In practice, we are actually interested in performing a 4 + (D — 4) split of the D dimensional 
space-time. This may be done as follows. The metric on a unit S D ~ 3 may always be written in 
terms of the line element on a unit S ^ 4 , denoted by diljj^4, as follows, 

hff ; (/.,-',/.,-' = d9 2 + sin 2 OdVlD-i , (2.11) 

where 9 is a polar- like coordinate, 9 £ [0, 7r]. Now we introduce four dimensional coordinates, 
:z+ = (a+, 9), fi — 0,1, 2, 3, and define a four dimensional metric 

g^dx^dx" = g^dx^dx" + f(x^)d9 2 , (2.12) 

as well as a new conformal factor 

A(x ft ) = sin 2 9g gg . (2.13) 
Then, the most general D dimensional metric compatible with SO(D — 2) isometry is, for D > 5 

ds 2 = g^dx^dx" + X(x^)dn D ^ . (2.14) 

Without specifying ([2TT2J) and ([2TT5J) . the geometry ([2~T4]) has only a manifest SO(D - 3) 
symmetry. We now perform a dimensional reduction on a (D — 4)-sphere. This yields, from 
the D dimensional vacuum Einstein equations, a set of 3 + 1 dimensional Einstein equations 
coupled to quasi-matter. If SO(D — 2) is the full isometry group, the quasi-matter terms do 
not contain independent degrees of freedom; rather, they may be completely determined by the 
3 + 1 dimensional geometry, via (|2.13|) . In this case, we could perform a dimensional reduction 
on a (D ~ 3)-sphcre, which has the full isometry group SO(D — 2). This would yield a 2+1 
dimensional system. The former method allows, however, the use of existing numerical codes, 
with small changes, which justifies our choice. 

The equations derived with dimensional reduction on a (D — 4)-sphere can be applied, of 
course, to describe also space-times in which the full isomctry group is SO(D — 3). This is the 
isomctry group of a class of BH collisions with impact parameter and with spin: the collisions in 
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which the two BHs always move on the same 2-plane and the only non trivial components of the 
spin 2-form arc on that same 2-plane - see Fig. [TJ With our framework we arc able, therefore, 
to describe not only head-on collisions of spinless BHs but also a class of collisions for spinning 
BHs with impact parameter. As follows from the discussion of (|2.9p . the ansatz (|2.14[) describes 
general space-times with SO{D — 3) isometry in D > 6. We remark that the models with D > 6 
are actually the most interesting for phenomenological studies of large extra dimensions models 
(see for instance (30^. 




FIG. 1: D dimensional representation, using coordinates (t , z), of two types of 

BH collisions: (left panel) head-on for spinless BHs, for which the isometry group is SO(D — 2); (right 
panel) non head-on, with motion on a single 2-plane, for BHs spinning in that same plane only, for which 
the isometry group is SO(D — 3). The figures make manifest the isometry group in both cases. 



2.2. Dimensional reduction on a (D — 4)-sphere and 3+1 split 

In the following we take (|2.14p as an ansatz, which has a manifest SO(D — 3) isometry. The 
D dimensional pure Einstein theory reduces then to a four dimensional theory of gravity coupled 
to a scalar field A(x M ). We remark that in this theory A and g^ are viewed as independent 
degrees of freedom; the relations (|2.12[) . (|2.13[) select a subset of the solution space. The solutions 
belonging to this subset have enhanced isometry SO(D — 2) and correspond to some of the 
physical processes we want to study (for instance, head-on collisions of spinless BHs). 

The D dimensional Einstein-Hilbcrt action reduces to 



1 

1u7g 4 



S = I d^x^—gX 



(D - - A -1 DA - — -^-X^d^Xd^X 



(2.15) 

where the D dimensional Newton's constant Gd is related to the four dimensional one G4 by 
the area of the unit D — 4 dimensional sphere: G4 = Gd/A s . Explicitly, the D dimensional 
Einstein's equations in vacuum yield the following system of four dimensional equations coupled 
to a scalar field: 



Rv» = ( V,AA - —d^XdvX ) , (2.16) 

VfyA = 2(D - 5) - ^— -dpWX. (2.17) 
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In these equations, all operators are covariant with respect to the four dimensional metric g^ u . 
The energy momentum tensor is 3 



T 



D-4 
16ttA 



l D — 5 
VjAA - — d^Xd^X - (D - 5).g M „ H j^g^d a Xd a X 



4A 



(2.18) 



With this four dimensional perspective, the usual 3 + 1 split of space-time [77|, l78j can be 
performed (see, e.g. [79|, H(|). For this purpose, we introduce the projection operator 7 P „ and 
the normal to the three dimensional hyper-surface E, n M (n^n^ = — 1), 

7/^f ~t f^pfov j 

as well as the lapse a and shift /3 M , 

<9 t = an + f3 , 

where t is the time coordinate. The four dimensional metric is then written in the form 



g^dx^dx" = -a 2 dt 2 + ^{dx 1 + ^dt^dx + f3 j dt) , i,j = 1,2,3. 



(2.19) 
(2.20) 
(2.21) 



As usual, we introduce the extrinsic curvature = —^C n jij, which gives the evolution 
equation for the 3-metric, 

(dt-Cp)'ri i =-2aK iJ . (2.22) 
The time evolution for Kij is given by 

{dt - Cp) = -Didja + a (WRij + KKij - 2K ik K k j) - a^tfjR^ , (2.23) 



where Di is the covariant derivative on the hyper-surface. The last term, 7 M ;7 I/ jR^v, vanishes 
for vacuum solutions. In the present case, it is given by the projection of equation (|2.16[) . 



D -4 
2A 



Using the formula 



and defining the variable 



we obtain 



D a D p X = -K^rfd^X + i^Yp^^X , 



K x = -\c n X = -^d ll X, 



VV/jV^A = DidjX - 2KijK x ■ 

Thus, (|2~2U)) becomes 

{d t - Cp) K^ = -Didja + a ^Rij + KK^ - 2K lk K k 



— a- 



D-A 
~2X~ 



DidjX - 2KijK x - -^piXdjX 



(2.24) 
(2.25) 
(2.26) 
(2.27) 

(2.28) 



3 We use the standard form of the Einstein equations G^v = 87rT M „ and choose geometrised units throughout. 
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To summarise, the evolution equations for the 3-metric and extrinsic curvature arc (|2.22[) and 

o. 

If the isometry group is SO(D — 3), the quasi- matter field A represents an independent degree 
of freedom, and we need to solve the evolution equations for A and K\. Even in the case of the 
larger isometry SO(D — 2), the evolution equations for A and K\ are useful as they enable us to 
test Eq. (|2.13[) and thus provide a check of the numerical evolution. The evolution equation for 
A is ([2^26)1 

(fit - Cp) A = -2aK x . (2.29) 
Eq. (|2.17[) provides an evolution equation for K\. The contraction of Eq. (|2.25[) with g a ", yields 

□ A = fWidjX - 2KK X - nW^A . (2.30) 



Noting that 



and 



we obtain 



C n K x = n^d^K x = -^V^cU - ^n v V^X , (2.31) 



»"V/ = -D"a , (2.32) 
a 



-rfn v V ud v \ = 2C n K x + -D v ad v X . (2.33) 

a 

Noticing also that D u ad v X = 7 U diadj A, we write 

□ A = --'P.iKX - 2KK\ + 2C n K x + -j^diadjX . (2.34) 
Moreover, from equation 

ZV = 7 %^A = d M X - 2n^K x , (2.35) 

we get 

d a \d a \ = f J d t \d ] \-AKl 1 (2.36) 
so that the evolution equation for Kx is 

- (d t - Cp) Kx = -^diXdja + (D-5) + KK X + ^r-^A'f - ^—^^d,Xd 3 X - ]-D k d k X . 
a la A 4A z 

(2.37) 

Equations (|2.29|) and (|2.37[) are the evolution equations for the quasi-matter degrees of freedom. 



2.3. BSSN formulation 



For numerical implementation, let us now write the evolution equations in the Baumgarte, 
Shapiro, Shibata and Nakamura (BSSN) formulation [6(1 [(J?}. Instead of evolving the variables 
jij and Kij , we introduce a conformal split of the physical 3-metric 7y as 

"fij = —% ■ (2.38) 
X 
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The conformal factor 

X = (det 7tf )" 1/3 , (2.39) 

is chosen such that det 7^ = 1 holds at all times. The extrinsic curvature is split into a conformal 
traceless part, A+j, and its trace, K, as 

Aj = x (K iS - ?fK) . (2.40) 

Moreover, we introduce the contracted conformal connection 

r = 7 Jfc f;. fc , (2.4i) 

where 

it r:; y x & kdjX + 5]kd%x ~ ™ kldi *) Tk = xfk + \^ ldix ' ^ 2 - 42 ) 

as an independent variable. In terms of the BSSN variables Xi Jij 1 A-ij ? T fc , the evolution equations 
are 

(ft - Cp) 7« = -2aA ij , (2.43a) 
(dt-£p)x=l<xxK, (2.43b) 
(d t - £p) K =[...}+ 4Tra(E + S) , (2.43c) 

(ft - £^) Ai S = [...]- 87ra ( X Sy ~ §7v) , (2.43d) 
(ft - £p) f* = [...] - lSTTa*" 1 j* , (2.43c) 

where [. . . ] denotes the standard right-hand side of the BSSN equations in the absence of source 
terms (see e.g. j!3|); the source terms are determined by 

E = n a n p T al3 , (2.44) 

ji = -ji a n T a p , (2.45) 

Stj = l a l l P 3 T aP , (2.46) 

S = f j Sij, (2.47) 

where the energy momentum tensor is given by Eq. (|2.18|) . A straightforward computation shows 
that 



4ir(E + S) 



-(D 5)\-i + ^A-V /2 7 tf A (x" 1/2 ftA) 



D-4 

~ 6 A- 2 xf 3 'ftAftA - \- x KK x - (D - h)X- 2 K 2 x , 



(2.48a) 



^ ^D- = ^A-'AftA + ^A- 1 (ftAftx + ftAftx - 7 W 7yftAft X ) - j%A" 2 ftAft -A 

- A- 1 /^^ - I^-A-^^V'A (x~ 1/2 ^a) + ^7«A- 2 X7 H ftAftA, 

(2.48b) 

16 ^ X _ * J> = 2 A : -"«•;.. A A - \- 2 K x ^d 3 \ - f^'^/A^ftA - ^-KX^djX, (2.48c) 



13 



where Di is the covariant derivative with respect to jij. 
Finally, the evolution equations for A and K\ are 



(8 t — Cp) A = -2aK x , (2.49a) 
h - C )K X =a\(D-5)+ [X-^BiXdjX - 4A~ 1 A^] 

(2.49b) 

KK\ \x* /2 l U Du 1 



As stated before, in the case of head-on collisions of spinless BHs the full symmetry of the 
D dimensional system we want to consider makes equations (|2.49[) redundant, by virtue of (|2.13[) . 
This allows to determine the quasi-matter degree of freedom in terms of the three dimensional 
spatial geometry, at each time slice. Indeed, we have only used an SO(D — 3) subgroup in the 
dimensional reduction we have performed. The extra symmetry manifests itself in the fact that 
jij possesses, at all times, (at least) one Killing vector field. If one chooses coordinates adapted 
to this Killing vector field, d/d6, the metric can then be written in the form (|2.12l) . and then the 
quasi-matter degree of freedom can be determined from the spatial geometry by (|2.13p . In the 
numerical implementation, one can either determine, at each time-step, the scalar field through 
(|2.13|) . or impose (|2.13[) only in the initial data, and then evolve the scalar field using Eq. (|2.49[) . 



3. INITIAL DATA 

Following the approach in (sTL |82| , we now derive the initial data of the evolution. 

3.1. D dimensional Hamiltonian and momentum constraints 

Let E be a (D - l)-dimcnsional space-like hyper-surface with induced metric and extrinsic 
curvature K ab in the D dimensional space-time. The space-time metric has the form 

ds 2 = g MN dx M dx N = -a 2 dt 2 + 7a6 (dx a + (3 a dt) (dx b + (3 b dt) , (3.1) 
where lower case latin indices take values a — 1, . ..,£) — 1. The constraint equations are 

R + K 2 - K ab K ab = , (3.2) 

D a (K ab - j ab K) = , (3.3) 

where R is the Ricci scalar of the hyper-surface S, K is the trace of the extrinsic curvature and 
D a is the covariant derivative with respect to j a b- 
We conformally decompose the spatial metric 

7a6 = ^ IT ^7a6, (3.4) 

which introduces the conformal factor tp, and split the extrinsic curvature in trace and trace- free 
parts, 

Kab = A ab + jy^jlab , (3.5) 

where "f ab A ab — 0. Define A ab = j ac j bd A ca <; define also the quantity 

> p+i 



A ab = ^— A ab , (3.6) 
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and lower its indices with jab, 

A ab = % c % d A cd = ^ 2 A ab . (3.7) 

Assuming that the "conformal metric" ^ a b is flat, which is a good approximation for the class 
of problems we want to study, we impose the "maximal slicing condition" K = 0. Then, the 
Hamiltonian and momentum constraints become 

V a i Q6 = 0, (3.8) 
M + 4 pl 3 2) ^-^A ab A ab = Q, (3.9) 

where V is the covariant derivative with respect to j a b and A is the flat space Laplace operator. 



3.2. Brill-Lindquist initial data and matching to four dimensions 

The simplest way to solve the constraints (j3.8p ~ p.9p is to require the extrinsic curvature to be 
zero 

K ab = 0. (3.10) 

This is sufficient to model the evolution of a single BH or even of N non-spinning, non-boosted 
BHs. The constraints reduce to a simple harmonic equation for the conformal factor, Atp = 0, 
which we solve in cylindrical coordinates {x a } = (z, p, 8, . . . ), where '. . . ' represent the coordi- 
nates on the (D — 4)-sphere, 

% b dx a dx b = dz 2 + dp 2 + p 2 (d9 2 + sin 2 6dQ, D -4) . (3.11) 

This choice of coordinates makes manifest the symmetries we want to impose. Observe that 9 
is a polar rather than an azimuthal coordinate, i.e. £ [0,7r]. Next, we introduce "incomplete" 
Cartesian coordinates as 

x = pcos9 , y = psinO 7 (3.12) 

where — oo < x < +oo and < y < +oo; we can then write the D dimensional initial data as 
p~TU]) together with 

% b dx a dx b = [dx 2 + dy 2 + dz 2 + y 2 dVL D ^] , (3.13) 

where ij) is a harmonic function on (|3 . 1 1 [) . 

If we compare the space-time metric (|3.1j) at the initial time slice, for which the spatial metric 
is given by (|3.4[) and (|3.13p . with the generic form that has an SO(D — 3) symmetry and is given 
by (|2.14p . (|2.2ip . we see that the initial data for the four dimensional variables are 

-..drdr - = ijjT^ [dx 2 + dy 2 + dz 2 ] , (3.14) 

and 

\ = y 2 ^nh>. (3.15) 

It remains to determine the initial conditions for Kij and K\. Using a set of D dimensional 
coordinates that make manifest the SO(D — 3) isometry, such as the one used in (|3 . 1 3|) . the 
vanishing of the extrinsic curvature is equivalent to 

#y = 0, (3.16) 
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whereas the vanishing of the components of K a b along the (D — 4)-sphcre implies that 

K x = 0. (3.17) 
Equations (|3.14|) f|3.17|l represent the Brill-Lindquist initial data in our framework. 

3.2.1. Evolution of a single black hole 

As one test of our framework we study the case of a single, non-spinning BH. Even though the 
space-time is static, the slicing evolves when using the puncture gauge. 

The solution for the conformal factor, which shall be used in the numerical tests to be presented 
below, is given by 

4 [x 2 + y 2 + {z - zbh) 2 } 

where the "puncture" [83[ is placed at x = y = and z = zbh- In this formulation, there is an 
interesting signature that the BH we wish to evolve is higher dimensional: the fall off of ip, which 
is that of a harmonic function in D — 1 spatial dimensions. Because the Tangherlini solution (69j 
may be expressed, in the same coordinate system as used in (|3.13|) . as 

/4R D - 3 - n D - 3 \ 2 / ,,D-3 \nh 

df = - { l R RD _ 3 + y 3 ) & + (l + ^zs) + dy* + dz 2 + yHn D . A ) , (3.19) 



where R = \/ x 2 + y 2 + z 2 , we conclude that the parameter /i appearing in the initial condition 
(|3.18p is the same which appears in this form of the Tangherlini solution. It is related to the 
ADM mass by 

D _ 3 = IftrAW (3 2Q) 

P (D~2)A SD ~ 2 { ' 

Note, however, that this form of the Tangherlini solution is not appropriate for a comparison 
with the numerical data. Indeed, the evolution does not, in general, preserve the conformally flat 
slicing of the initial condition, which is the slicing used in this form of the Tangherlini solution. 
We shall return to this issue in Section I4TT1 



3.2.2. Head-on collision of black holes 

As another test of our formulation, and in particular of the numerical code's long term stability, 
we also evolve a head-on collision of non-spinning non-boosted BHs. In this case, the initial data 
for the conformal factor arc given by 

u D - 3 ll d - 3 

^ = 1 + — TD-3U2 + — nrru2 ■ ( 3 - 21 ) 

4 [x 2 +y 2 + (z- z A ) 2 ] ( 3)/ 4 [x 2 +y 2 + (z- z B ) 2 ] ( }/ 



This conformal factor is used in Section POI 
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THE NUMERICAL TREATMENT 



Our numerical simulations have been performed by adapting the Lean code [6g, initially 
designed for 3+1 vacuum space-times. The Lean code is based on the Cactus computational 
toolkit [84]]. It employs the BSSN formulation of the Einstein equations [66,|67|, uses the moving 



puncture method 0, H| , the Carpet package for mesh refinement [85l. l86f. the spectral solver 
described in [gjt for 3+1 initial data and Thornburg's AHFinderDirect |88l.l89j. Details about 
Lean may be found in [68j . Here we focus on the numerical issues generated by the quasi-matter 
terms arising from the dimensional reduction by isometry. 

We expect that the quasi-matter field A has a y 2 fall off as y — > 0, that is, on the xz plane. 
This leads to divisions by zero on the right-hand side of the BSSN evolution equations, cf. (|2.48[) . 
Since we expect all variables to remain regular on the xz plane, all divisions by y need to be 
cancelled by a corresponding fall off behaviour of the numerators. At y = 0, however, in order 
to implement this behaviour numerically, we need to isolate the irregular terms and evaluate 
expressions such as 

lim-, (4.1) 

y^O y 

where / is some example function which behaves like y n with n > 1 near the xz plane. It is 
necessary, for this purpose, to formulate the equations in terms of variables which are manifestly 
regular at y = 0. We also prefer to apply a conformal re-scaling of A and use the evolution 
variable 

C^4 A - ( 4 - 2 ) 

As in (|2.49p , in order to obtain a first order evolution system in time, we introduce an auxiliary 
variable (see Appendix [A]) : 

= "5^a (ft " £/0(Cl/ 2 ) = ~ (dtC - P m d m ( + \td m r - . (4.3) 

2ay z 2a \ 3 y J 

The third term on the right-hand side arises from the fact that ( is not a scalar, but a scalar 
density of weight —2/3. The inclusion of this term might not be necessary for a stable numerical 
implementation. For consistency with the rest of the BSSN variables, however, we decide to keep 
this form of Kq. 

The quasi-matter terms (|2.48[) , the quasi-matter evolution equations (|2.49[) and the constraints 
are recast in terms of £ and Kq in Appendix [21 In particular, we notice that 

K X = V -K C + \V^K. (4.4) 

X 3 x 

A detailed analysis of the equations in terms of the variables £ and Kq shows how all terms with 
an explicit dependence on l/y n , n > 1 may be treated for numerical implementation. This is 
discussed in Appendix iBl 



4.1. Numerical results in D — 5 



Wc first address the question of longevity of our simulations in D = 5. It is also of interest 
in this context to test the code's capability to successfully merge a BH binary. For this purpose 
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we have evolved a head-on collision starting from rest. The initial conditions are those from 
Section IXO with 

A = Ml = y , (4-5) 
z A = -z B = 3.185 n, (4.6) 

and we use the grid setup (cf. Sec. II E of Ref. J68|) 

{(512, 256, 128, 64, 32, 16, 8) x (2, 1), h = 1/32} , 

in units of fx. The gauge variables a and f3 l are evolved according to the modified moving puncture 
conditions (|A5[) and (|A6[) with parameters t]k = r\Kc = 1-5 an d ?y = 0.75. We employ fourth 
order discretization in space and time and impose a floor value [2[ for the variable x = 10~ 4 . 




FIG. 2: The BSSN variable \ e h panel) and the quasi-matter momentum K( (right panel) are shown 
along the axis of collision for a head-on collision at times t = 0, 5, 20, 40 and 256 fi. Note that Kq = 
at t = 0. 



In Fig. [5] we show the conformal factor \ an d the momentum Kq along the axis of collision 
at various times. At early times, the evolution is dominated by the adjustment of the gauge 
(cf. the solid and short-dashed curves). The two holes next start approaching each other (long- 
dashed and dotted curves) and eventually merge and settle down into a single stationary hole 
(dash-dotted curves). We have not observed any signs of instability and decided to stop the 
simulation at t = 256 /i. It is reassuring to notice that the framework can handle the merger 
in as robust a fashion as has been demonstrated by various numerical groups for BH binaries in 
3+1 dimensions. 

We have also used the head-on collision to test the relation between the scalar field A and the 
3 + 1 metric discussed in Sec. 11.31 for the case that SO(D — 2) is the full isometry group. We 
have verified for this purpose that Eq . 12 . 1 31 remains satisfied to within a relative error of 10~ 3 in 
the immediate vicinity of the puncture and at most 10 ~ 5 everywhere else. 

In order to further test our numerical framework, we have performed simulations of a single 
BH, using the initial data described in Section f3 . 2 . 1 1 and the grid setup 

{(512, 256, 128, 64, 32, 16, 8, 4, 2) x (), h} , 

in units of \i with resolutions h c = 1/32 and hf = 1/48. In Fig. [3] we show the Hamiltonian 
constraint and the y component of the momentum constraint at evolution time t = 28/i. By 
this time there is hardly any more gauge dynamics going on. One can see that there is some 
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noise, but the overall convergence is acceptable. For the Hamiltonian constraint the convergence 
is essentially 4th order and for the momentum constraint it decreases slightly towards 2nd or 
3rd order in patches. From experience in 3+1 dimensional numerical relativity this is perfectly 
acceptable, especially given the fact that prolongation in time is second-order accurate. 




0.01 0.1 1 10 100 0.01 0.1 1 10 100 

y/ji y/ji 
(a) Hamiltonian constraint (b) y-component of the momentum constraint 



FIG. 3: Constraints at time t = 28/i, for the evolution of a single Tangherlini BH in five dimensions. 



A different test of our numerical code was performed in order to compare the analytical 
Tangherlini solution with our numerical results. The challenge to do this comparison, at the 
level of the line element, is to write the well known analytical solution in the same coordinate 
system in which the numerical evolution is occurring. One way around this problem is to fix the 
numerical gauge as to match a known coordinate system for the analytic solution. Following (64j 
we fixed the gauge parameters to be 

a = l, f3 l = 0, i = 1,2,3; (4.7) 

this corresponds to geodesic slicing. The D dimensional Tangherlini solution may be expressed 
in a coordinate system of type (|3.ip with a = l 1 f3 a = Q,a = l,...,D — 1. This coordinate system 
may be achieved by setting a congruence of in-falling radial time-like geodesies, each geodesic 
starting from rest at radial coordinate ro, with ro spanning the interval [/i, +oo[, and using their 
proper time r and ro as coordinates (instead of the standard t, r Schwarzschild-like coordinates). 
A detailed construction of the Tangherlini solution in these coordinates is given in Appendix [C] 
The line element becomes 

(ro(R) 2 + (t^)) 2 >V dR 2 / ( V X 

ds 2 = -dr 2 + A \ T^~w+ -[77m) t2 W' ^ 

where r (i?) is given by Eq. (|C5[) . 

The numerical evolution in this gauge is naturally doomed. Geodesies hit the physical singular- 
ity at finite proper time. Thus, this slicing is inappropriate for a long term numerical evolution. 
As long as the evolution does not break down, however, there is perfect control over the slicing, 
and hence the numerical and analytical evolution can be compared with ease. This is shown in 
Fig. 01 where we have plotted one metric component ^ xx along the x axis (left) and £/x (right), 
for various values of r using both the analytical solution and numerical data. The agreement is 
excellent for ^ xx and good for C/x- The latter shows some deviations very close to the puncture, 
but we believe that it is not a problem for two reasons: (i) the agreement improves for higher 
resolution; (ii) the mismatch does not propagate outside of the horizon. 
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(a) ^xx along the ai-axis 



(b) C/X = ^/V 2 along the j/-axis 



FIG. 4: Numerical values versus analytical plot (solid lines) for various values of r, for the single 
Tangherlini BH in five dimensions. The horizontal axes are labelled in units of /i. 

It is easy to interpret the behaviour observed for ^ xx . The geodesic that starts from r — tq (in 
Schwarzschild-like coordinates) hits the physical singularity of the Tangherlini solution within 
proper time r = r\j [i. Moreover, this happens at 



R 



(4.9) 



The earliest time at which the slicing hits the singularity is r = which happens at R = /i/2. 
On the x-axis R = x and indeed one sees in Fig.2]that j xx diverges at x = /j,/2. The divergence 
then extends to both larger and smaller values of x, as expected from 



4.2. Preliminary numerical results in D = 6 

A quick glance at the evolution equations (|A7aP and (|A7b[) of the scalar field £ as well as the 
source terms (|A8a[) - ([A8c[) indicates that D = 5 may be a special case. In all these expressions 
there exist terms which manifestly vanish for D = 5. In contrast, there exist no terms which 
manifestly vanish for any dimension D > 6. The purpose of this Section is to extend the test of 
our framework to a case which involves all source terms. 

We have indeed noticed one fundamental difference between simulations in D = 5 and those 
using D > 6. Whereas we have been able to obtain stable simulations of single BHs lasting 
hundreds of \i for the former case by modifying the moving puncture gauge conditions, we have 
not yet succeeded in doing so for D > 6. While the lifetime of the simulations in D > 6 shows 
a dependence on the exact nature of lapse and shift, all simulations developed instabilities on 
a timescale of about 10 \i. Resolving this issue requires an extensive study involving a large 
number of experiments with gauge conditions, constraint damping and possibly other aspects 
of the formulation. Such a study is beyond the scope of this work and deferred to a future 
publication. The results presented in this Section still provide valuable information. Most 
importantly, they demonstrate the internal consistency of the code for D > 6 and thus minimise 
the possibility of a simple error in the implementation. Furthermore they exhibit clearly that 
our framework and in particular our regularisation of the variables as discussed in Appendix [B] 
is in principle suitable for simulations in arbitrary dimensions. 

Wc first consider the convergence of the constraints analogous to the results displayed in Fig. [3] 
for D = 5. Compared to those simulations, the only change we have applied in D = 6 is to set 
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the gauge parameters to t\k = ?7if c = f] — 2. This choice enables us to evolve single BHs to 
about 10 fi when instabilities cause the runs to abort. In Fig. [5] we show the Hamiltonian and 




y/fi 

(a) Hamiltonian constraint (b) y-component of the momentum constraint 

FIG. 5: Constraints at time t = 8^i, for the evolution of a single Tangherlini BH in six dimensions. 



the y-component of the momentum constraint at t = 8 [i along the y-axis. As for D — 5, the 
high resolution result is amplified by a factor 1.5 4 expected for fourth order convergence [80l |. 
While the convergence appears to be closer to second order in some patches of the momentum 
constraint, the results are clearly compatible with the numerical discretization. 

For the second test, we compare the numerical evolution of a single D = 6 Tangherlini BH 
with the analytic solution, using geodesic slicing. This comparison is more difficult in the present 
case than in D = 5, because the line element analogous to (|4.8p cannot be obtained in a simple 




FIG. 6: Numerical values versus the semi-analytic expression of % 
for the single Tangherlini BH in six dimensions. 



(cf. Appendix [C]) along the z-axis 



analytic form. In Appendix [Cl we demonstrate how a semi-analytic solution can be obtained for 
the metric. In Fig. [6] we compare this expression with the three dimensional numerical values at 
times r = 0.5 fj,, 0.7 fj, and 0.72 [i. The agreement is excellent and demonstrates that our code 
works well at least up to the point where instabilities set in. As mentioned above, resolving these 
stability problems will be of the highest priority in future extensions of our work. 
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5. FINAL REMARKS 

In this paper we present a framework that allows the generalisation of the present generation of 
3+1 numerical codes to evolve, with relatively minor modifications, space-times with SO(D — 2) 
symmetry in 5 dimensions and SO(D — 3) symmetry in D > 6 dimensions. The key idea is 
a dimensional reduction of the problem along the lines of Geroch's t 70j procedure that recasts 
the D dimensional Einstein vacuum equations in the form of the standard four dimensional 
equations plus some quasi-matter source terms. The resulting equations can be transformed 
straightforwardly into the BSSN formulation that has proved remarkably successful in numerical 
evolutions of BH configurations in 3+1 space-times. We have isolated several issues related to 
the regularisation of the variables used in our formulation and demonstrated how all difficulties 
related to the coordinate singularity arising out of the use of a "radius-like" coordinate can be 
successfully addressed in a numerical implementation. We have further illustrated how initial 
data for single, non-spinning BHs as well as BH binaries with vanishing initial extrinsic curvature 
can be adapted straightforwardly to the formulation presented in this paper. More generally, the 
class of problems that may be studied with our framework includes head-on collisions in D > 5 
and a subset of BH collisions with impact parameter and spin in D > 6. 

As might be expected, stable evolutions of such space-times require some modifications of the 
underlying methods of the so-called moving puncture technique, especially with regard to the 
gauge conditions used therein. We have successfully modified the slicing condition via incorpo- 
ration of the canonical momentum of the quasi-matter field in order to obtain long-term stable 
simulations in D = 5 dimensions. Unfortunately, these modifications do not appear sufficient to 
provide long-term stability for arbitrary values of the dimensionality D. We will address this 
important issue in the form of a systematic study in future work. 

Wc have tested our framework by adapting the Lean code and performed a variety of single BH 
space-times. Most importantly, we have demonstrated the internal consistency of our numerical 
framework in D = 5 and 6 dimensions by showing convergence of the Hamiltonian and momentum 
constraints as well as comparing numerical results with (semi-)analytic expressions for a single 
Tanghcrlini BH in geodesic slicing. We have further shown for D = 5 that the head-on collision 
of a BH binary successfully merges into a single hole which settles down into a stationary state 
and can be evolved numerically for long times, hundreds of fi in the present example. 

A complete study of such BH binary evolutions requires the implementation of gravitational 
wave extraction in arbitrary dimensions as well as the generalisation of apparent horizon diag- 
nostics beyond D = 4. Both are currently being implemented in the Lean code and will be 
discussed in detail in future work. 

In spite of several open questions, we believe that our formalism will open up a vast range 
of uncharted territory in BH physics for contemporary numerical relativity. The list of possible 
applications and extensions of our framework is too large to be included here, and we merely 
mention strong hyperbolicity studies of the BSSN formulation with sources and systematic in- 
vestigation of BH binary dynamics in D dimensions. These studies are under way and will be 
reported elsewhere. 
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Appendix A: Implementing regular variables 

The numerical evolution faces a problem at the symmetry axis, given the quasi-matter terms in 
(I2.48|) and the initial data discussed in Sec. [3] The "incomplete" Cartesian coordinate y vanishes 
at the symmetry axis, cf. (|3~T2| . Then, from (|3~T5j) . A vanishes at the axis (except, possibly, at 
the puncture). Inspection of equations (|2.48[) and (|2.49[) immediately reveals various divisions 
by A, leading to numerical problems. 

From previous experience with polar and spherical coordinates in simpler models involving, 
for example, neutron stars (cf. [90, |91|) we know that it is better to avoid the use of singular 
variables such as A. We should use, instead, regular functions. In our case, since A behaves as 
y 2 near the axis, this is simply achieved by introducing a variable C via (|4.2p . The evolution of 
£ is formulated in terms of a first order in time system of equations. For this purpose we have 
introduced in Eq. (|4.3p the variable Kq. We remark that if, instead, we employ the standard 
definition for the momentum associated with £, i.e. 

K C = ~^(dt- C P )C = ~ (d t C - p m 8 m ( + 2 -Qd m p m ^j , (Al) 

we face problems in the numerical evolution for vanishing lapse. This may be seen as follows. 
From (|4.4|) 

K x =V-K i + \^K + ^. (A2) 

X 3 x ax 

Acting on both sides of this equation with the derivative operator 

d = dt - P m dm , 



results in the expression 

"2- 



d K c = ^0 K X + 4!—K c + -aKK c + -(—K + -(aK 2 



(A3) 

- ^d K + - — - -K c d m p m + C — 5 d ° a ■ 

3 a \ y J 6 a y ya z 

This is an evolution equation for K^. To obtain it explicitly one uses (|2.37[) to express doK\, 
together with 

d K = -D m d m a + a (.A mn A mn + ^K 2 ) + 4™(£ + S) . (A4) 
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Moreover we need gauge conditions. Throughout this work we use the following coordinate 
choices 

d a a = -2a(r] K K + Vk ( K c ) , (A5) 
d P = - Ty/f . (A6) 

Note the extra term involving Kq in the slicing condition compared with standard moving punc- 
ture gauge in 3+1 dimensions and the additional freedom we have introduced in the form of the 
parameters t\k and r\K, ■ 

The problems in the case of a collapsed lapse become clear if we consider the final two terms 
in (|A3|) . These terms do not change when BSSN variables are introduced and diverge for the 
modified moving puncture gauge conditions (|A5[) and (|A6[) as the lapse a — > 0. We have solved 
this problem by expressing our equations in terms of the variable Kq (|4.3[) . instead of Kq. 

In BSSN variables, the evolution equation for £ and Kq (which replace the quasi-matter evo- 
lution equations (|2.49[1 ) become 

d t C = -2aK ( + f3 m d m ( - l(d m p m + 2C— , (A7a) 

3 y 

d t K c = p m d m K c - \K Q d m r + 2^K C - \(d K - ^ m d m a - \l mn {d m a){ X dnC - (d nX ) 
6 y 6 y A 

(5 - D)^(^y - 1) + (4 - ™cU + ^^-i ym d m x 
y y 2 y 

4 Q 4 4 X 

+ {D- 6)^ + ^-^KK ( + ^j^CK 2 + \i m% ((D m d nX - X D m d n () +xCy 

(A7b) 

These equations have no manifest problems as a — > 0. 
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In terms of the regular variables, £ and K^, the quasi- matter terms (|2.48l) read 



Att{E + S) 
D-A 



D-A 



Y ^vvr _ i 2D - 7 T y 

(D-^^ ^-i mn {d m O{d nX )-x- 



2C 



4 C 2 

KKt 1 o W n (y 

3 y \C 



C 

D — l 



8ttx (Sij - ^7ijS) _ 1 



-1 



4' 1 



, {dmX)(dnX) 
X 



d m Q - d m x 



K c K 



1 jV™ 
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9 m X 



■2\ 



D i d jX + ^-7ij7 m "5„x dmX - T<9mC - lij d m x 



TF 



Kr 1 



7 J 



167rji _ 2 
£-4 ~ y 



1 



1 



i 



i 



Finally, the constraints are now given by 



(A8a) 



(A8b) 



(A8c) 



H = R + -K 2 - r in l kl A m kA n i - IGirE , 



Mi = T 



DA --A ^ 

2 X 



where we also need to express E in terms of our fundamental variables. It is given by 



(D - 3 )4f m S r „C -(D- 2)-^ m d mX + ^-^l mn {d m Q{d n O 

yQ y 4 c 



£)-4 



(A9) 
(A10) 



D-2 



D + 3_, 



2C A\ C 



Kt 2D-A„K e 



-K- 



£) + 1 o y ~ T v y y yv C - 1 

-^-K 2 + ±r m Dmd n ( - i mn D m dnX ~ 2x— + {D - 5)^ 7 \ 

9 C y C y 2 



(All) 



Appendix B: Analysis of troublesome terms at y = 

The right-hand sides of Eqs. (|A7|) - (|A10[) contain various terms which cannot be evaluated 
directly at y — because they involve explicit division by y. Although these terms are regular 
by virtue of a corresponding behaviour of the numerators, they need to be explicitly evaluated 
in the numerical implementation. In this Appendix we outline how the regularity of these terms 
can be implemented in a simple and efficient manner. For convenience we use a special notation: 
late latin indices i, j, ... run from 1 to 3, covering x, y and z, but early latin indices a, b, 
take values 1 and 3 but not 2, i.e. they cover x and z but not y. 
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We begin this discussion by describing a simple manipulation which underlies most of our 
rcgularisation procedure. Consider for this purpose a function h which is linear in y near y = 0, 
i.e. its Taylor expansion is given by h(y) = h\y + 0(y 2 ). From this relation we directly obtain 

lim - = hi = d v h . (Bl) 

y^O y 

This trading of divisions by y for partial derivatives extends to higher orders in a straightforward 
manner and will be used throughout the following discussion. 

Next, we consider the right-hand sides of Eqs. (|A7[) - (|A10[) and summarise the potentially 
troublesome terms as follows 

RV TV 

-, (B2) 

y v 

~yrn 

—d m f, (B3) 

y 

:LJ ^> B4 
2T 

~ (Si y ^-l ym A m ?J , (B5) 

i(«C + ^5 i C-2Cff J .) . (B6) 

Here / stands for either of the scalars or densities C, \ an d a - 

Regularity of the terms (|B2[) immediately follows from the symmetry condition of the y- 
componcnt of a vector 

P y {-y) = -P y (y)- (B7) 
We can therefore use the idea illustrated in Eq. (|B1[) and obtain 

lim = d y (3 v , (B8) 

y^O y 

and likewise for T y jy. The terms (|B3[) are treated in a similar manner because the derivative of 
a scalar (density) behaves like a vector on our Cartesian grid. We thus obtain 



lim ^d m f = (d y j ya )(daf) + J yy dyd y f . (B9) 



Regularity of the expression (|B4[) is not immediately obvious but can be shown to follow 
directly from the requirement that there should be no conical singularity at y = 0. Specifically, 
this condition implies that ^ vv ( = 1 + 0(y 2 ), so that 

'jyy(-i\ _ 1 



lh a ( ) = o ({Wr™ + "f yv W) ■ (bio) 



The discussion of the term (|B5[) requires us to distinguish between the cases i = a ^ y and 
i = y. The former straightforwardly results in 

lim (-T^A ma ) = -A ba d y 7 yb - i yy d v A ya . (Bll) 



2G 



For the case i = y, we first note that the limit y — >• implies j yy = l/% y + 0(y 2 ), so that the 



(B12) 



condition (|B10I) . i.e. no conical singularities, can be written as 

lim(C-%,)=0(y 2 ). 

y->0 

Next we take the time derivative of this expression and obtain after some manipulation 



0(y 2 ) = hm ft(C - i vy ) = -2a( - ^ m A my 

V-*0 \ Q 



0(y 2 ) 



and, consequently, 



lim 

y-s-0 



1 K, 



y 



= o. 



(B13) 



(B14) 



Finally, we consider the term (|B6|) . Expansion of the Christoffel symbol, repeated use of the 
method illustrated in Eq. (|Blj) and the condition for avoiding a conical singularity enable us to 
regularise this term for all combinations of the free indices i and j. We thus obtain 



lim 

y-s-o 

lim 



lim 



-2-r 



2d y d y ( - Cl vv d y dy% y - ad v */ vc )(2d v % c - d c % y ) 
0, 

-C,^ VV {d y dalby + dyd b % a ~ dydylab) 

- ({d y j yc )(d a j bc + d b % c - d c % b ) . 



(B15) 
(B16) 

(B17) 



We conclude this discussion with a method to express derivatives of the inverse metric in terms of 
derivatives of the metric. For this purpose we use the condition that det jij = 1 by construction 
and explicitly invert the metric components as for example in 



T V = Ixzlyz - Ixylzz 



(B18) 

A straightforward calculation gives us the derivatives of the inverse metric components as follows 

dyT V = lxzdy%z ~ lzzdy%y + O ^ ) , (B19) 



dyl VZ = Ixzdylxy ~ Ixxdylyz + 0(y 2 ) , 

dy^ VV = J ZZ dy% X + Jxxdyjzz ~ 2^ xz 3y^ XZ , 

dydyl VV = lz Z d v d y % x + % x d y d y % z - 2% z d y d y % z + 0(y 2 ) 



(B20) 
(B21) 
(B22) 



The benefit in using these expressions is purely numerical: we do not need to store the inverse 
metric in grid functions which reduces the memory requirements of the simulations. 



Appendix C: Geodesic slicing 



In standard Schwarzschild-like coordinates, the Tangherlini metric reads 



ds 2 



-f{r)dt 2 



dr 2 

m 



r 2 dn 



D-2 



m = 1 - (~ 



D-3 



(CI) 
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For a radially in- falling massive particle, starting from rest at r — ro, the energy per unit mass 
is \J /(?'())■ The geodesic equation may then be written as 



dt = yTfo) 



(C2) 



In four and five dimensions these equations have simple solutions. In five dimensions the solutions 
are 



t= VfW)r+^ln 



+ Vf(ro)r 2 /n 



^1 

r 0, 



Then, performing a coordinate transformation (t, r) — > (r, ro) the line element becomes 



(C3) 



ds 2 



-dr 2 



r l + I 7T, 



dr 2 



rg/(ro) 



(C4) 



This coordinate system encodes a space-time slicing with zero shift and constant (unit) lapse 
{i.e. of type (|3.1[l with a = l,/3 a = 0) for all times. To compare it with a numerical evolution 
we must have the initial data for the spatial metric written in a conformally flat form. Taking 
the initial hyper-surface to be r — we see that this is achieved by a coordinate transformation 
ro — > R with 



dR 
R 



drn 



\/f(ro)r 



r (R)=R[ 1 



4i? 2 J 



(C5) 



This actually coincides with the standard coordinate transformation from Schwarzschild to 
isotropic coordinates in five dimensions. The line element finally reads (|4.8[) . At the initial 
hyper-surface r = 0, 



ds r=o = 



MR) 
R 



(dR 2 + R 2 dQ 3 ) = 



r (y/p 2 +z 2 ) 
^p 2 + z 2 



(dz 2 + dp 2 + P 2 d0 2 + p 2 sin 2 9d£l x ) , 



where we have used the metric on the 3-sphere in the form 

dfl 3 =d8 + sin 2 8(d8 2 + sin 2 0dfi x ) , 
and performed the coordinate transformation (R, 6) — > (p, z) defined as 

p = Rsin 9 , z = R cos 6 . 

Using (|3~T2]) we get 



ds 



T=0 



ro(V x 2 + V 2 + z ' 2 ) 
\/ x 2 + y 2 + z 2 



(dx 2 + dy 2 + dz 2 + y 2 dVti) 



(C6) 
(C7) 

(C8) 
(C9) 



Thus the coordinate transformation from the spherical coordinates used in (|4.8[) . (R, 9, 8), to the 
"incomplete" Cartesian coordinates used in the numerical evolution (x, y, z) is 



x = R sin 9 cos 9 , y = R sin 8 sin 9 , z = R cos 9 , 



(CIO) 



2s 



which resembles the usual coordinate transformation from spherical polar coordinates to Carte- 
sian coordinates in R 3 ; but note that 9 and 9 are both polar angles with range [0,7r], which is 
the manifestation of the Cartesian coordinates "incompleteness" . 

The coordinate change (|C10|) brings the five dimensional Tangherlini metric in geodesic slicing 
to a conformally flat form at r = 0. This matches the initial data for the numerical evolution. 
One may ask, however, if the coordinate transformation evolves, in order to compare the analytic 
form with the numerical evolution. This cannot be the case, since the existence of r-dependent 
terms in the coordinate transformation would imply a drift away from geodesic slicing. We are 
thus guaranteed that the coordinate transformation (|C10[) is valid for all values of r. Then, 
we can predict the value of the metric components that should be obtained from the numerical 
evolution; say j xx should be, at time r 

j xx {r,x,y,z)- R2 + R i( x 2 + y 2 ) + ( x 2 +2/2)2 . lyu) 

where R 2 = x 2 + y 2 + z 2 and 3_r_r(t, R),ggg(r, R),ggg(r, R) are readily obtained from (|4.8[) with 
(|C7|) and (|C10[) . The result for j xx along the x-axis is plotted in Fig. H] for various values of r. 

For D > 6 the situation is more involved because equations (|C2|) can no longer be integrated 
straightforwardly, but require a numerical treatment. First one notices that the coordinate 
transformation (t, r) — > (r, tq), with initial conditions t(r = 0) = and r(r = 0) = tq, brings the 
D dimensional Tangherlini metric to the form 

ds 2 = -dr 2 + ( ^ZllAY^L + r 2 (Tj ro)< m D _ 2 . (C12) 



dr ) f(r ) 

Then, from the initial conditions, it follows that the coordinate transformation to isotropic 
coordinates at r = is 

dR D ^ r„(R) = Z(l + jLY / 



r a{ R) = -ll + J— . (C13) 

R v/M^-o A 4 V 4i ?V 

Writing the metric on the (D — 2)-sphere as in (|C7|) (replacing dCli —> dVLn^^), one concludes 



that the transformation to "incomplete" Cartesian coordinates is still (jClOp . Thus (|C11I) is still 
valid, which reduces to, along the x-axis (R = x): 

/ nnl / s r (x) 2 / dr(T,r )\ 2 
j xx {t,x,0,0) = g RR {T,x) = —[ — . (C14) 

X \ Or /r =ro(x) 

This expression is valid for any D. For D = 6, ro(x) is explicitly given by (|C13|) . The derivative 



in (|C14[) has to be computed numerically. The result for j xx is plotted, for various values of r, 
in Fig. |B1 
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